Edwards thermodynamics of the jamming transition for frictionless packings: 
ergodicity test and role of angoricity and compactivity 
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This paper illustrates how the tools of equilibrium statistical mechanics can help to explain a 
far-from-equilibrium problem: the jamming transition in frictionless granular materials. Edwards 
ideas consist of proposing a statistical ensemble of volume and stress fluctuations through the ther- 
modynamic notion of entropy, compactivity, X, and angoricity, A (two temperature- like variables). 
We find that Edwards thermodynamics is able to describe the jamming transition (J-point). Using 
the ensemble formalism we elucidate the following: (i) We test the combined volume-stress ensemble 
by comparing the statistical properties of jammed configurations obtained by dynamics with those 
averaged over the ensemble of minima in the potential energy landscape as a test of ergodicity. 
Agreement between both methods supports the idea of "thermalization" at a given angoricity and 
compactivity. (ii) A microcanonical ensemble analysis supports the idea of maximum entropy prin- 
ciple for grains. (Hi) The intensive variables describe the approach to jamming through a series of 
scaling relations as A — > + and X — Y CP . Due to the force-volume coupling, the jamming transition 
can be probed thermodynamically by a "jamming temperature" Tj comprised of contributions from 
A and X. (iv) The thermodynamic framework reveals the order of the jamming phase transition by 
showing the absence of critical fluctuations at jamming in observables like pressure and volume, (v) 
Finally, we elaborate on a comparison with relevant studies showing a breakdown of equiprobability 
of microstates. 
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The application of concepts from equilibrium statisti- 
cal mechanics to out of equilibrium systems has a long 
history of describing diverse systems ranging from glasses 
to granular materials [IH1]- For dissipative jammed 
systems — particulate grains or droplets — the key con- 
cept proposed by Edwards is to replace the energy en- 
semble describing conservative systems by the volume 
ensemble Q. However, this approach alone is not able 
to describe the jamming point (J-point) for deformable 
particles like emulsions and droplets [4j-|7| , whose geomet- 
ric configurations arc influenced by the applied external 
stress. Therefore, the volume ensemble requires augmen- 
tation by the ensemble of stresses [§4Tl| . Just as volume 
fluctuations can be described by compactivity, the stress 
fluctuations give rise to an angoricity, another analogue 
of temperature in equilibrium systems. 

In the past 20 years since the publication of Edwards 
work there has been many attempts to understand and 
test the foundations of the thermodynamics of powders 
and grains. Three approaches are relevant to the present 
study: 

1. Experimental studies of reversibility. — Start- 
ing with the experiments of Chicago which were 
reproduced by other groups fl2T - fl5l |. a well-defined 
experimental protocol has been introduced to 
achieve reversible states in granular matter. These 
experiments indicate that systematically shaken 
granular materials show reversible behavior and 
therefore are amenable to a statistical mechan- 
ics approach, despite the frictional and dissipative 
character of the material. These results are com- 



plemented by direct measurements of compactiv- 
ity and effective temperatures in granular media 
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2. Numerical test of ergodicity. — Numerical sim- 
ulations compare the ensemble average of observ- 
ables with those obtained from direct dynamical 
measures in granular matter and glasses. These 
studies [19l - l24j find general agreement between 
both measures and, together with the experimental 
studies of reversibility fl2l - fl5j , suggest that ergod- 
icity might work in granular media. 

3. Numerical and experimental studies of 
equiprobability of jammed states. — Exhaus- 
tive searches of all jammed states are conducted 
in small systems to test the equiprobability of 
jammed states, as a foundation of the microcanon- 
ical ensemble of grains. Numerical simulations and 
experiments indicate that jammed states are not 
equiprobable [25l - [27| . These results suggest that a 
hidden extra variable [28[ might be needed to de- 
scribe jammed granular matter in contrast with the 
work in 1 and 2. 

The current situation can be summarized as follow- 
ing: When directly tested or exploited in practical ap- 
plications, Edwards ensemble seems to work well. These 
include studies where ensemble and dynamical measure- 
ments are directly compared, and recent applications 
of the formalism to predict random close packing of 
monodisperse spherical particles [29l [30| . polydisperse 
systems [311 ], and two [32| and high dimensional systems 
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33]. However, a direct count of microstates reveals prob- 
lems at the foundation of the framework manifested in 
the breakdown of the flat average assumption in the mi- 
crocanonical ensemble (25l - [28| . 



In this paper we investigate the Edwards ensemble 
of granular matter focusing on describing the jamming 
transition [j]-|7[. A short version of this study has been 
recently published in [34[ ■ We employ a strategy that 
mixes the approaches 2 and 3 above. We first perform 
an exhaustive search of all jammed configurations in the 
Potential Energy Lands cape (P EL) of small frictionless 
systems in the spirit of |25l - t2a |. We then use this infor- 
mation to perform a direct test of ergodicity in the spirit 
of [UHl- Ovrr results indicate: (i) The dynamical and 
ensemble measurements of presure, coordination number, 
volume, and distribution of forces agree well, supporting 
ergodicity. A microcanonical ensemble analysis supports 
also a maximum entropy principle for grains, (ii) In- 
tensive variables like angoricity, A, and compactivity, X, 
describe the approach to jamming through a series of 
scaling relations. Due to the force-volume coupling, the 
jamming transition can be probed thermodynamically by 
a "jamming temperature" Tj comprised of contributions 
from A and X. (Hi) These intensive variables elucidate 
the thermodynamic order of the jamming phase transi- 
tion by showing the absence of critical fluctuations above 
jamming in static observables like pressure and volume. 
That is, the jamming transition is not critical and there 
is no critical correlation length arising from a thermody- 
namic n-point correlation function. We discuss other pos- 
sible correlation lengths, (iv) Surprisingly, we reproduce 
the results of [HI regarding the failure of equiprobability 
of microstates while obtaining the correct dynamics mea- 
surements as in [lQi-24]. We then offer a possible solu- 
tion to this conundrum to elucidate why the microstates 
seems to be not equiprobable while the ensemble averages 
produce the correct results. 



The paper is organized as follows. Section Q] discusses 
the Edwards thermodynamics of the jamming transition. 
Section |TT] describes the ensemble calculations in the Po- 
tential Energy Landscape formalism. Section Mil de- 
scribes the Hertzian system to be studied. Section IIVI 
describes the ensemble measurements to be compared 
with the MD measures of Section [V] Section PVTl explains 
how to calculate A from the data. The ergodicity test is 
made in Section I VIII Section IVIIII describes the calcula- 
tion in the microcanonical ensemble where the principle 
of maximum entropy is verified and the coupled jamming 
temperature is obtained. Section HXl compares our results 
with those of O' Hern et al. [25[ and Section [X] summa- 
rizes the work. Appendix |A1 includes "de yapa" a study of 
coordination number fluctuations in the Edwards theory 
for random close packings of hard spheres. 



I. EDWARDS THERMODYNAMICS AND THE 
JAMMING TRANSITION 

The process typically referred to as the jamming tran- 
sition occurs at a critical volume fraction <j) c where the 
granular system compresses into a mechanically stable 
configuration in response to the application of an external 
strain |2|, |J] . The application of a subsequent external 
pressure with the concomitant particle rearrangements 
and compression results in a set of configurations char- 
acterized by the system volume V — NV g /<fi ((f> is the 
volume fraction of N particles of volume V g ) and applied 
external stress or pressure p (for simplicity we assume 
isotropic states). 

It has been long argued whether the jamming transi- 
tion is a first-order transition at the discontinuity in the 
average coordination number, Z, or a second-order tran- 
sition with the power-law scaling of the system's pressure 
as the system approaches jamming with 4>-4> c ^0+0- 
0, HH . Previous work [ll|, [3(1 H3] has proposed to explain 
the jamming transition by a field theory in the pressure 
ensemble. Here, we use the idea of "thermalization" of 
an ensemble of mechanically stable granular materials at 
a given volume and pressure to study the jamming tran- 
sition from a thermodynamic viewpoint. 

For a fixed number of grains, there exist many jammed 
states [IE HH confined by the external pressure p in a 
volume V. In an effort to describe the nature of this 
nonequilibrium system from a statistical mechanics per- 
spective, a statistical ensemble [1, H3, El was introduced 
for jammed matter. In the canonical ensemble of pres- 
sure and volume, the probability of a state is given by 
exp[-W(dS/dV) - T(dS/dT)], where S is the entropy 
of the system, W is the volume function measuring the 
volume of the system as a function of the particle coor- 
dinates and r = pV is the boundary stress (or internal 
virial) [36] of the system. Just as dE/dS = T is the 
temperature in equilibrium system, the temperature-like 
variables in jammed systems are the compactivity Q 

X = dV/dS, (1) 

and the angoricity 0, 

A = dT/dS. (2) 

In a recent series of papers [29M33| the compactiv- 
ity was used to describe frictional and frictionless hard 
spheres in the volume ensemble. Here, we test the valid- 
ity of the statistical approach in the combined pressure- 
volume ensemble to describe deformable, frictionless par- 
ticles, such as emulsion systems jammed under osmotic 
pressure near the jamming transition 38, 39]. 

In general, if the density of states ^(r, </>) in the space 
of jammed configurations (defined as the probability of 
finding a jammed state at a given (r, cf>) at A = oo) is 
known, then calculations of macroscopic observables, like 
pressure p and average coordination number Z as a func- 
tion of 0, can be performed by the canonical ensemble 
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average [36|, |37| at a given volume: 
1 



pg(r,0) e" ai dr, (3) 



and 



(Z(a,<f>)) ena = ^J o Z g(T, 



e - aT dr, 



(4) 



where the canonical partition function is 

poo 

z= / 5 (r» e - Qr dr, (5) 

Jo 

and the density of states is normalized as <?(r, <p)dT = 
1. The inverse angoricity is defined as 



a=l/A = dS/dT. 



(6) 



At the jamming transition the system reaches isostatic 
equilibrium, such that the stresses are exactly balanced 
in the resulting configuration, and there exists a unique 
solution to the interparticle force equations satisfying me- 
chanical equilibrium. It is well known that observables 
present power-law scaling 0-0|: 



(p)dyn ~ (0 - <j>cT , 



(7) 

(z) dyn -z c ~ {4>-4>c)\ (8) 

where a — 3/2 and b = 1/2 for Hertzian spheres and 
Z c = 6 is the coordination number at the frictionless 
isostatic point (J-point) The average (--^dyn in- 

dicates that these quantities are obtained by averaging 
over packings generated dynamically in either simula- 
tions or experiments as opposed to the ensemble average 
over configurations (• • • ) en s of Eqs. ©-([4]). Comparing 
the ensemble calculations, Eqs. ©-(H]), with the direct 
dynamical measurements, Eqs. ([?])— (|8]l. provides a basic 
test of the ergodic hypothesis for the statistical ensemble. 

Our approach is the following: We first perform an 
exhaustive enumeration of configurations to calculate 
g(T, 0) and obtain (p(a, 4>)}cns a s a function of a for a 
given cj) using Eq. ([3]). Then, we obtain the angoric- 
ity by comparing the pressure in the ensemble average 
with the one obtained following the dynamical evolution 
with Molecular Dynamics (MD) simulations. By setting 
{p(a, (j)))ens = (p)dyn, we obtain the angoricity as a func- 
tion of <f>. By virtue of obtaining a(cf>), all the other ob- 
servables can be calculated in the ensemble formulation. 
The ultimate test of ergodicity is realized by comparing 
the remaining ensemble observables with the correspond- 
ing direct dynamical measures. 



II. POTENTIAL ENERGY LANDSCAPE 
APPROACH: ENSEMBLE CALCULATIONS 

A. Features of the Potential Energy Landscape 

An appealing approach for understanding out-of- 
equilibrium systems is to study the properties of the 



system's "potential energy landscape" (PEL) [4J|, de- 
scribed by the 3A-coordinates of all particles in the 
multi-dimensional configuration space, or landscape, of 
the potential energy of the system (N is the number 
of particles). Characterizing such potential energy land- 
scapes has become an important approach to study the 
behaviour of out-of-equilibrium systems. For example, 
this approach has provided important new insights into 
the origin of the unusual properties of supercooled liq- 
uids, such as the distinction between "strong" and "frag- 
ile" liquids [42| . 

In frictionless granular matter, the potential energy is 
well-defined and each jammed configuration corresponds 
to one local minimum in the PEL. For small systems 
(N iS 14), it is possible to find all the minima with cur- 
rent computational power [25| . For somewhat larger sys- 
tems N « 30, it is possible to obtain a representative 
ensemble, without exhaustively sampling all the states. 
Based on these stationary points, we test the combined 
volume-stress ensemble. The following work is only valid 
for frictionless systems where the potential energy of in- 
teraction is well defined. Frictional grains are path de- 
pendent due to Coulomb friction between particles and 
therefore not amenable to a PEL study since there is no 
well defined energy of interaction. 

The formalism introduced by Goldstein 43j consists 
of partitioning the potential energy surface into a set of 
basins as illustrated in Fig. [1] The dynamics on the po- 
tential energy surface can be separated into two types: 
the vibrational motion inside each basin and the transi- 
tional motion between the local minima. Stillinger and 
coworkers [44[ developed the method of inherent struc- 
ture to characterize the PEL. In this method, a local 
minimum in the PEL is located by following the steepest- 
descent pathway from any point surrounding the mini- 
mum. The inherent structure formalism simplifies the 
energy landscape into local minima and ignores the vi- 
brational motion around them. The dynamics between 
the inherent structures is introduced with the transition 
states identified with the saddle points in the PEL. The 
transition states are stationary points like the local min- 
ima but they have at least one maximum eigendirection. 



B. Finding Stationary States 

For the simplest system of N structureless frictionless 
particles possessing no internal orientational and vibra- 
tional degrees of freedom, the potential energy function of 
this N-body system is E(ri, . . . , rjv), where the vectors 
7"j comprise position coordinates. As mentioned above, 
the most interesting points of a potential energy surface 
are the stationary points, where the gradient vanishes. 
Here we explain how to locate these stationary points. 
The algorithm follows well established methods in com- 
putational chemistry |4lj . The procedure is analogous 
to finding the inherent structures [45[ of glassy systems. 
The algorithm employed, LBFGS algorithm, is also sim- 
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FIG. 1: A model two-dimensional potential energy surface. 
The energy landscape is divided into basins of attraction, 
where the minima are the jammed states connected by path- 
ways through saddle points. States A and B are typical pack- 
ing configurations of 30 particles (in blue) with their periodic 
boundary systems. 

ilar to the conjugate gradient method employed by O' 
Hern [25|, [2(1], differing in the fact that it does not 
require the calculation of the Hessian matrix at every 
time step. We make the source code in C++ available 
at http://www.jamlab.org and free to use together with 
all the packings generated in this study. The algorithm 
has been used in the short version of this article (34J 
and in a study of the PEL in Lennard- Jones glasses to 
reconstruct a network of stationary states and apply a 
percolation picture of the glass transition (4(| . 

C. General Method — Newton- Raphson Method 

Consider the Taylor expansion of the potential energy, 
E, around a general point in configuration space, r, 

E(r + h) = E{r)+g T h+-h T Hh + 0(h 3 ), (9) 

where g is the gradient, gi = diE, H is the Hessian ma- 
trix, H{j = didjE, and h is a small step vector that gives 
the displacement away from r. 

By Eq. ©, the calculation of energy difference for a 
given step h from the initial point r is complicated. By 
selecting the eigenvectors of the Hessian matrix e a as our 
local coordinates, we can simplify the Taylor expansion 
of Eq. © as: 

AE = E(r + h)-E(r)^J2(9 a h a + ^h 2 a ), (10) 

a 

where g = Y. a 9o<e a , h = ^ a /i a e Q , He a = X a e a , and 
X a is the eigenvalue of the Hessian matrix for component 
a. 
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FIG. 2: A schematic energy change curve for one component 
with X a > 0. We can select the downhill step as h a — — j^- 
to obtain a maximum energy change. The uphill step can not 
be too large since the Taylor expansion will not be accurate 
enough for the calculation. Here, the uphill step is chosen as 

From Eq. (TTOj) . it is easy to see that the total change 
of energy could simply be the sum of the changes in each 
directions. This may help us to raise the energy in some 
directions and reduce the energy at others, and finally 
reach a stationary point. The length of each step compo- 
nents can be selected as the maximum change of energy: 

h a = Sa-r^, (11) 

as shown in FigjS] The sign S a — ±1 in this formula 
depends on the choice of uphill or downhill direction. In 
fact, for X a > 0, it is possible to choose another step for 
the uphill case, since AE a increases as \h a \, but for large 
steps, the Taylor expansion Eq. ([9]) may breakdown. 
Therefore, it is important to control the step length. For 
X a < 0, we reach the opposite conclusion. 

The stationary points can be separated into local min- 
ima and saddle points. Based on the eigenvalues of the 
Hessian matrix for the stationary point, the local minima 
are ordered as: 

< Ai < A 2 ■ • • < A 3 at, (12) 

and for a saddle point of order a: 

Xi < ■■■ < X a <0 < A Q+1 < ••• < X 3N . (13) 

Generally, this algorithm searches for the nearest sta- 
tionary point on the surface by following the opposite 
(A a > 0) and along (A Q < 0) the various gradient direc- 
tions. 



D. Finding local minima — LBFGS algorithm 

It is much easier to locate local minima than saddle 
points because, for the first, we only need to search down- 
hill in every direction. At present one of the most effi- 
cient methods to search the local minima for large system 
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is Nocedal's limited memory Broyden-Fletcher-Goldfarb- 
Shanno algorithm (LBFGS) [ME!- The LBFGS algo- 
rithm constructs an approximate inverse Hessian matrix 
from the gradients (first derivatives) which are calculated 
from previous points. Since it is only necessary to cal- 
culate the gradients at each searching step, the LBFGS 
algorithm increases the computational speed of the algo- 
rithm enormously. 

In the Newton-Raphson method discussed above, 
the Hessian matrix of second derivatives is needed 
to be evaluated directly. Instead, the Hessian 
matrix used in LBFGS method is approximated 
using updates specified by gradient evaluations. 
The LBFGS algorithm code can be obtained from 
|http: // www. netlib . org / opt / index. html| Here we present 
a brief explanation of the algorithm 

From an initial random point ro and an approximate 
Hessian matrix Hq (in practice, Hq can be initialized 
with Hq = /), the following steps are repeated until r 
converges to the local minimum. 

• Obtain a direction h k by solving: H k h k — 
-VE(r k ). 

• Perform a line search to find an acceptable step 
size 7^ in the direction found in the first step, then 
update r k+ i = r k + j k h k . 

• Set s k = a k h k . 

• Set y k = VE(r k+1 ) - VE(r k ). 

T 

• Set the new Hessian, H k +\ = H k + Vk J )k — 

H k s k (H k s k ) T 
s k H k s k 




Local Minimum A Local Minimum B 



FIG. 3: A two dimensional 31 particle system in a circular 
boundary. Three different configurations in this system are 
generated with different algorithms. The LBFGS method is 
applied to locate minima A and B. For saddle C which con- 
nects A and B, the eigenvector following method is used. 



where A Q are the eigenvalues of the Hessian matrix and 
g a are the components of the gradient in the diagonal 
base (h a is set to for the directions where X a = 0). 
The sign S a = ±1 is chosen by the order of the saddle 
point. For a saddle point of order n, the algorithm will 
set S a = — 1 for 1 < a < n and S a — 1 for a > n. 

When g a —> 0, the step size h a converges to the 
Newton-Raphson step as Eq. (|11[) : 

h a = S a ^ + 0(g 2 a ). (15) 



F. An Example 



E. Finding saddles — Eigenvector following method 

In the present study we do not make use of the saddle 
points. However, other studies using network theory to 
represent the PEL necessitate the links between minima 
through the saddle points [46|. For completeness, below 
we explain how to search for saddles. A particular pow- 
erful method for locating saddle points is the eigenvector 
following method [41| . 

The eigenvector-following method, developed by Cer- 
jan, Miller and others [U I4814521 ] . consists of locating a 
saddle point from a local minimum. At each searching 
step towards a saddle point with a order, the directions 
are separated into two types: a uphill directions to maxi- 
mization and 3N — a downhill directions to minimization. 

We follow the implementation of the eigenvector- 
following method by Grigera [4!|. We give a general 
description: at each searching step, a step size h is cal- 
culated by the diagonalized Hessian matrix [49|, [5l|, HI] : 



S a - 



\ a \(i + VT+M/% 



(14) 



We generate a two dimensional soft-ball system in cir- 
cular boundary, which contains 31 particles of equal ra- 
dius, to illustrate the method of finding stationary and 
saddle points in the PEL. The interaction between par- 
ticles (also the interaction between particles and wall) 
follows the Hertzian law [6|: 



V(n,rj) = e\n - rj\'- 



(16) 



Here, e is the interaction strength between particles i and 
j. The volume fraction is 4> — 0.80, which is closed to 
the jamming transition of 2d hard disks. 

We first generate a random configuration, which is 
the initial point of the search of the minima. With the 
LBFGS method, we search the local minimum A nearby 
this initial point. After the minimum A is obtained, we 
apply the eigenvector following method to walk from the 
point A on the potential energy surface to locate the tran- 
sition state C (here the transition state is a first order 
saddle). Finally, the minimum B is located by applying 
LBFGS method again. Figure [3] shows configurations of 
two local minima (marked as red) and the transition state 
(marked as blue) between them. 
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FIG. 4: The pathway from minimum A to minimum B, pass- 
ing by the saddle C, the x-coordinate is the distance from sad- 
dle C, the y-coordinate is the potential energy of the packing. 

The pathway from minimum A to minimum B, passing 
by transition state C, is shown in Fig. HJ The pathway 
distance is the Euclidean distance, 

d = yj{f- r)(r' - r) = - r it<x f, (17) 

y i,ot 

where i = 1,2,3, a = 1---3N, r' is the coordinate of 
configuration passing along the searching method and r 
is the coordinate of saddle C. 

The dynamics from minimum to minimum can be rep- 
resented as a walk on a network whose nodes corre- 
spond to the minima and where edges link those min- 
ima which are directly connected by a transition state. 
The work of Doye |53f provides an illustration of such a 
landscape network for a LJ energy surface. To charac- 
terize the topology of the landscape network, Doye [53[ 
study small Lennard- Jones clusters to locate nearly all 
the minima and transition states on the potential en- 
ergy landscape. The inherent structure network of such 
a system has a scale-free and small-world properties. In 
a companion study [46j we repeated the main results as 
Doye studied. The numbers of minima and transition 
states are expected to increase roughly as N min ~ e aN 
and N s t ~ Ne aN respectively, where TV is the number 
of atoms in the cluster. Therefore, the largest network 
that we are able to consider is for a 14-atom cluster for 
which we have located 4158 minima and 90 738 transition 
states in agreement with the results of Doye. In the next 
Section we apply the above formalism to find the station- 
ary states for a 3d granular system of Hertz spheres in a 
periodic boundary. 

III. SYSTEM INFORMATION. HERTZIAN 
SYSTEM OF SPHERES 

Next we calculate the density of jammed states g(T,<p) 
in the framework of the PEL formulation for a system of 



Hertz spheres. In the case of frictionless jammed systems, 
the mechanically stable configurations arc defined as the 
local minima of the PEL [{J [26| . 

The systems used for both, ensemble generation and 
molecular dynamic simulation, are the same. They are 
composed of 30 spherical particles in a periodic boundary 
box. The particles have same radius R — 5/im and inter- 
act via a Hertz normal repulsive force without friction. 
The normal force interaction is defined as @, HH, HH : 

F n = H k n R^ 2 (6r) s , (18) 

where 6 = 3/2 is the Hertz exponent, Sr — (1/2) [2R — 
\xi — a?2 1] > is the normal overlap between the spheres 
and k n = 4G/(1 — v) is defined in terms of the shear 
modulus G and the Poisson's ratio v of the material from 
which the grains are made. We use typical values for 
glass: G = 29 GPa and v = 0.2 and the density of the 
particles, p = 2 x 10 3 kg/m 3 @, [H[. The interparticle 
potential energy is 

E = ljTT Rl/2{Sr)S+1 - (19) 

The Hertz potential is chosen for its general applica- 
bility to granular materials. The results are expected to 
be independent of the form of the potential. Below, we 
apply the LBFGS algorithms [H, S3 to find the local 
minima of the PES (zero-order saddles). 

IV. ENSEMBLE GENERATION 

In this section, we first explain the method to obtain 
geometrically distinct minima in the PEL to calculate 
the density of states. Then we show that the density 
of the states, g(T,(f>), does not change significantly after 
sufficient searching time for the configurations. 

In principle, if all local minima corresponding to the 
mechanically stable configurations of the PEL are ob- 
tained, the density of states g(T,<fi) can be calculated. 
Such an exhaustive enumeration of all the jammed states 
requires that N not be too large due to computational 
limits. On the other hand, in order to obtain a precise 
average pressure in the MD simulation, (p)d yn , N cannot 
be too small such that boundary effects are minimized. 
Considering these constraints, we choose a 30 particle 
system. 

In order to enumerate the jammed states at a given vol- 
ume fraction <j>, we start by generating initial unjammed 
packings (not mechanically stable) performing a Monte 
Carlo (MC) simulation at a high, fixed temperature. The 
MC part of the method applied to the initial packings 
assumes a flat exploration of the whole PEL. Every MC 
unjammed configuration is in the basin of attraction of 
a jammed state which is defined as a local minimum in 
the PES with a positive definite Hessian matrix, that 
is a zero-order saddle. In order to find such a minimum, 
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we apply the LBFGS algorithm provided by Nocedal and 
Liu [43] explained above. The PEL for each fixed <f> likely 
includes millions of geometrically distinct minima by our 
simulation results. Therefore, an exhaustive search of 
configurations is computationally long; for a system of 
30 particles it is impossible to find all the configurations 
with the current available computational power. How- 
ever, we notice that it not crucial to find all the states, 
but rather a sufficiently accurate density of states. There- 
fore, we check that the number of found configurations 
has saturated after sufficient trials and that the density 
of states g(r,<f>) has converged to a final shape under a 
prescribed approximation. 

It is also important to determine if the local minima 
are distinct. Usually, the eigenvalues of the Hessian ma- 
trix at each local minimum can be used to distinguish 
these mechanically stable packings. Here, we follow this 
idea to compare minima to filter the symmetric pack- 
ings. However, instead of calculating the eigenvalues of 
each packing, which is time consuming, we calculate a 
function of the distance between any two particles in the 
packing to improve search efficiency (for the LBFGS al- 
gorithm, we do not need to calculate the Hessian matrix). 
For each packing, we assign the function Qi for each par- 
ticle: 



Qi 



E 

1<3<N, jjH 



9 

7rr ■ 
tan 2 (-^), 



3L 2 



(20) 



where is the distance between particles i and j, L is 
the system size and N = 30. We list the Qi for each 
packing from minimum to maximum {Qi}(l < i < N). 
Since Qi is a higher order nonlinear function, we can 
assume that two packings are the same if they have the 
same list. The tolerance is defined as: 



T = 



'Er 



<i<N 



(Qi 



N 2 



(21) 



where Qi and Q i are the corresponding values from the 
lists of two packings. 

Figure [5] shows the distributions of the tolerance T 
for packings at different volume fractions. This figure 
suggests that two packings can be considered the same if 
T < 10 -1 , which defines the noise level. 

From Fig. |51 we see that after one week of search- 
ing, g(T, (p) does not change significantly, since the initial 
packings are generated by a completely random protocol. 
We also calculate the probability of finding new mechani- 
cally stable states for different searching days, defined as 
N new (i) /N t otai(i) , where N ncw (i) is the number of new 
configurations found on the i-th day and iVtotai(i) is the 
total number of configurations found in i days. From Fig. 
El we see that, after one week searching, the probability of 
finding new configurations at different volume fractions 
seems to have converged in the linear plot. Figure [TJd 
shows a detail of the actual number of new configurations 
found and g(T,<fi) versus searching time in days suggest- 
ing convergence. However, the log-log plot of the inset in 




logT 



FIG. 5: The distribution of the tolerance T between any two 
packings at the given <j>. From the graph, the value of T for 
which any two different packings are considered to be same 
is chosen to be 10 _1 , which is above the noise threshold and 
below the distribution of T. 



10' 



10' 
10' 



10-' , 



[^^^^^ ( a ) 


Day 1 

Day 2 

Day 3 

Day 4 

Day 5 
Day 6 
Day 7 
Day 8 
Day 9 

Day 10 








Day 11 

Day 12 

Day 13 

Day 14 

Day 15 



FIG. 6: Log-log plot of the distribution of g(T, (j>) for 15 
searching days (a) at <j> — 0.609, (b) at <j> = 0.614, (c) at 
4> = 0.625. Different color in (a), (b), (c) corresponds to the 
different day. We find that after 15 days the distributions 
have converged. 



Fig. [T^i indicates that the algorithm is still searching for 
new configurations; the power-law relation in the inset 
suggesting a neverending story. However, the main ques- 
tion is whether the observables have converged. A further 
test of convergence is obtained below in Fig. [T4l where the 
value of the inverse angoricity is measured as a function 
of the searching time in days. This plot suggests that 
enough ensemble packings have been obtained to cap- 
ture the features of g(T, </>) that give rise to the correct 
observables. We conclude that we have obtained an ac- 
curate enough density of states for this particular system 
size. Regarding system size dependence, the presented 
results are still N dependent, although they started to 
converge for N ~ 35 and above, Fig. [5] More accurate 
calculations for large values of N remain computation- 
ally impossible, but in our treatment the exact choice of 
N is not as important as the consistency of the results 
between ensemble and MD, for a given N value. 
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FIG. 7: (a) The probability to find new configurations as a 
function of searching time, (b) Linear plot of the density of 
states as a function of searching time. Different colors indicate 
different days according to the inset. Inset shows the actual 
number of new configurations. 



Figure |H] shows g(T, <j)) versus T for different volume 
fractions. 



V. MD CALCULATIONS 

In order to analyze numerical results, we perform MD 
simulations to obtain Zd yn and </>dyn, which are herein 
considered real dynamics. The algorithm is described in 
detail in [29|, HH, IHl • Here, a general description is given: 
A gas of non-interacting particles at an initial volume 
fraction is generated in a periodically repeated cubic box. 
Then, an extremely slow isotropic compression is applied 
to the system. The compression rate is To = 5.9tg , 
where the time is in units of to — -Ry p/G. After obtain- 
ing a state for which the pressure p is a slightly higher 
than the prefixed pressure we choose, the compression is 
stopped and the system is allowed to relax to mechani- 
cal equilibrium following Newton's equations. Then the 
system is compressed and relaxed repeatedly until the 
system can be mechanically stable at the predetermined 
pressure. To obtain the statical average of Zd yn and </>dyn, 
we repeat the simulation to get enough packing samples 
having statistically independent random initial particle 
positions. Here, 250 independent packings are obtained 
for each fixed pressure (see Fig. [TU]). (j) — (</>)dyn and 




15 20 25 30 35 40 
System size 

FIG. 8: Dependence of the results on the system size. The 
average value of p converges as early as N ~ 25 particles. The 
distribution g(F, (f>) (inset) has not fully converged yet but its 
shape has converged after N = 35 and the first moment does 
not change as indicated by the average p. 
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FIG. 9: The density of states g(T, cj>) as a function of internal 
virial T for different volume fraction, cj>, ranging from 0.610 to 
0.670. The inset shows the logarithmic distribution of g(T, <j>). 
At low volume fraction (</> < 0.625), the distributions are 
sharp and the tails of the distributions are exponential. At 
high volume fraction (<j> > 0.640), the distributions are much 
broader and the tails are Gaussian. 



(Z)dyn are flat averages of these 250 packings by 

Sl<i<250 & 



and 



(^dyn 



( Z )dyn 



250 



zCl<i<250 Z i 

250 ' 



(22) 



(23) 



From previous studies, it has been observed the pres- 
sure p vanishes as power-law of <p when approaching the 
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FIG. 10: The cyan O is ( a ) 0dyn and (b) Za yn for every 
single packing obtained with MD and the blue O is (<^)dyn 
and (^)dyn average over the single packings for the system 
which are shown in the text of the paper. 



jamming transition as seen in Eq. (0 [H, [|| . We obtain 

(Fig. ED 

<p)dy„=P0 (cWc) 1 - 68 , (24) 

where <p c = 0.6077 is the volume fraction corresponding 
to the isostatic point J @, @| following Eq. ([5} and p a = 
10.8MPa. This critical value (j> c and the exponent, a — 
1.65, are slightly different from the values obtained for 
larger systems (a — 5) 0, 0]. However, our purpose is 
to use the same system in the dynamical calculation and 
the exact enumeration for a proper comparison. 



VI. ANGORICITY CALCULATION 

Since we obtain g(T, <p) and (p)dyn for each volume 
fraction 0, we can calculate the inverse angoricity a by 
Eq. ([3]). The pressure (p(a, <p))ens for a given ^ is a 
function depending on a as: 



(p(a, 0)) ons = 



f o °° ra (r,0)e"« r dr _ J2pe 



aT 



g(T,<f>)e-*r d r E« 



-aT 



(25) 



Figure [12] shows the result of the numerical integration 
of Eq. (|2"5")) for a particular <fi = 0.614 as a function of 




FIG. 11: Scaling of pressure. The blue O shows the power- 
law relation for (p)d yn vs (0)dyn — <t>c for the 30-particle system. 
Here, the pressure (p)d yn are average values obtained by 250 
independent MD simulations. The red O is the pressure used 
to obtain the inverse angoricity a predicted by Eq. (|24|) . 
The relatively small system size results in large fluctuations 
of the observables. In order to predict a precise relation for 
the system (N = 30), sufficient independent samples of the 
packings are generated to calculate the precise average for 
observables. We prepare 250 independent packings for each <f> 
to get enough statistical samples to obtain (p)dyn and {Z)d yn 
by statistical average. The inset shows a semi-log plot. 
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FIG. 12: The numerical integration of Eq. (25JI for <j> = 0.614 
is shown as the pink curve. We input the (p)dyn (pink O i n 
the plot) and obtain the corresponding inverse angoricity a. 



a using the numerically obtained g(T, <f) from Fig. |H1 
To obtain the value of a for this 4>, we input the corre- 
sponding measure of the pressure obtained dynamically 
(p(4>))d yn and obtain the value of a as schematically de- 
picted in Fig. [T^] The same procedure is followed for 
every <j> (see Fig. [T3|) and the dependence a(<fi) is ob- 
tained. 

We also check the inverse angoricity ol(4>) using g(T,<p) 
for different searching days, to ensure the accuracy and 
convergence to the proper value. From Fig. Q3J we can 
see that, after 10 days searching, a(4>) is stable due to the 
fact that the density of state, g(T, <j>), does not change 
significantly. 

For each <j> we use g(T,(f>) to calculate (p(a)) ons by 
Eq. ©. Then, we obtain a(<j>) by setting {p(a,4>)) cns — 
(p)dyn for every cf). The resulting equation of state a{(f>) is 
plotted in Fig. [15] and shows that the angoricity follows 
a power-law, near <p c , of the form: 



A oc 



(26) 
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FIG. 13: Calculation of a for several volume fractions <f> as 
explained in detail in Fig. 1121 
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FIG. 14: Calculation of inverse angoricity a as a function of 
searching time. 
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FIG. 15: (a) Inverse angoricity a as a function of 4>-<t>c- We 
find a power-law relation for system's volume fraction <j> near 
4>c- The solid line has a slope of -2.5. (b) The angoricity 
A(= 1/a) vs (f>-(f>c- To find A accurately for system's volume 
fraction (j> much larger than <f> c , becomes difficult due to the 
large fluctuations and finite size effects. In principle, we ex- 
pect that the plateau of A for large volume fraction <f> might 
be related to the finite size of the sample. Indeed it is very 
difficult to estimate a since it falls in the plateau in Fig. 1131 



with 7 = 2.5. The result is consistent with 7 = S + 1.0, 
suggesting that A oc T oc F n r. For volume fraction much 
larger than 4> c , the system's input pressure (p(</>))dyn 
reaches the plateau at low a of the function (p(a, 4>)} en s 
(see Fig. I13j) and the corresponding a(<p) becomes much 
smaller (the angoricity A((f>) becomes much larger), lead- 
ing to large errors in the value of A as 4> becomes 
large. This might explain the plateau found in A when 
(cj) — 4> c ) > 2 x 10~ 2 as shown in Fig. [15] 

Angoricity is a measure of the number of ways the 
stress can be distributed in a given volume. Since the 
stresses have a unique solution for a given configuration 
at the isostatic point, (f> c , the corresponding angoricity 
vanishes. At higher pressure, the system is determined by 
multiple degrees of freedom satisfying mechanical equi- 
librium, leading to a higher stress temperature, A. The 
angoricity can also be viewed as a scale of stability for the 
system at different volume fractions. Systems jammed at 
larger volume fractions require higher angoricity (higher 
driving force) to rearrange. 



VII. TEST OF ERGODICITY 

In principle, using the inverse angoricity, a, from Eq. 
(pMf we can calculate any macroscopic statistical observ- 
able (fi)em at a given volume by performing the ensemble 
average |37j : 

1 I" 00 

(s(0)) e „s = zJ a B 5(r, 4>) e- aT dr. (27) 

We test the ergodic hypothesis in the Edwards's ensemble 
by comparing Eq. (|27[) with the corresponding value ob- 
tained with MD simulations averaged over (250) sample 
packings, -B,, generated dynamically: 

250 

W))dy„ = ^]>>" ( 28 ) 
i— 1 

The comparison is realized by measuring the average 
coordination number, (Z), the average force and the dis- 
tribution of interparticle forces. We calculate (Z) cus by 
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Eq. (@J and (Z) dyn as in Eq. (|28[) . Using cv(</>) for each 
volume fraction, we calculate (Z) eas by: 



TO)>« 



Jo°° Zg(T, , 



,-aT 



dr _ J2ze 



-ctT 



The average force (F) cns is given by: 



(29) 



TO))ens = 



/ o oo Fg(r,0)e-» r dr _EF e 

/ o oo 5 (r,0)e-«r d r E« 



(30) 



where F is the average force for each ensemble packing. 
Finally, the force distribution P cns (F/F) is given by: 



Pens(F/F) 



f™P(F/F)g(T,0)e- ar dT _ £P(F/F)e 



0(r, 



-ar dr 



(31) 

Equations (f29 |) - (|3"Tj) are then compared with the dynam- 
ical measures for a test of ergodicity in Figs. [TBI and HPTl 
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(Z) ens = (Z)dyn- The average inter-particle force F for 
a jammed packing is proportional to the pressure of the 
packing. We calculate (F) cns and (-F)dyn and find that 
they coincide very closely (see Fig. UTah The full dis- 
tribution of inter-particle forces for jammed systems is 
also an important observable which has been extensively 
studied in previous works [5(| ■ The force distribu- 
tion is calculated in the ensemble P cns (F/F) by averaging 
the force distribution for every configuration in the PES. 
Figure ITTb shows the distribution functions. The peak 
of the distribution shown in Fig. [TTb indicates that the 
systems are jammed d, [H, Besides the exact shape 
of the distribution, the similarity between the ensemble 
and the dynamical calculations shown in Fig. [TTb is sig- 
nificant. The study of (Z), (F) and P(F/F) reveals that 
ar the statistical ensemble can predict the macroscopic ob- 
'servables obtained in MD. We conclude that the idea of 
"thcrmalization" at an angoricity is able to describe the 
jamming system very well. 
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FIG. 16: Test of ergodicity. (a) The blue O is the av- 
erage coordination number {Z)d yn obtained by 250 indepen- 
dent MD simulations. The red O i s the coordination number 
{Z) ena calculated by the ensemble for different volume frac- 
tions. Agreement between both measures supports the con- 
cept of ergodicity in the system. (b)The same as (a) but 
in a log-log plot. The blue O shows the power-law rela- 
tions for {Z)d yn -Z c vs (4>)dyn ~4>c for 30-particle system with 
4> c = 0.6077 and Z c = 5.82. 

Figure [Trlh and ITBb show that the two independent 
estimations of the coordination number agree very well: 
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FIG. 17: Test of ergodicity. (a) Comparison of (i ? )dyn and 
{F) cns for different volume fractions, (b) The comparison of 
selected distribution of force P dyn (F/F) and P en s(F/F) for 
different volume fractions. 

The MD simulations performed so far are at a prede- 
termined pressure p. For this case there is no difference 
between the force distribution P(F/F) and P(F/(F)) @. 
On the other hand, a MD simulation at a given fixed vol- 
ume fraction (j), gives rise to different distributions. For 
each system with fixed (j), the packings can have various 
pressure. This suggests that the force distribution for 
each packing scaled by the average force over all pack- 
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ings, P(F/(F)), should be different from the force dis- 
tribution scaled bythe average force of that particular 
packing P{F/F) [25|]. We now proceed to investigate a 
constant volume MD, vMD simulation. 



Pcns{F/ (F) cns ) in the ensemble average: 




FIG. 18: (a) The distribution of force P vM d(F/ (F) vM d)- (b) 
The distribution of force P C ns(F/(F) cns ). (c) and (d) The 
comparison of selected P(F/ (F)) between vMD and ensemble 
predicted by angoricity. 



P ens (F/(F) 

I. 



f °°P(F/(F) cns )g(rA)t 



dr 



,-or 



dr 



-, (32) 



where (F) cns is the overall average F of the ensemble. 

From Fig. [T5b . we find the same tendency as obtained 
in MD simulation. Furthermore, we check the distribu- 
tion of force P(F/(F)) for our vMD system (see Fig. 
TO]) . We see that P(F/(F)) for different volume fraction 
<j> collapses very well similarly to those obtained from the 
predetermined pressure system. This result suggests that 
P(F/ (F)) is a global quantity that can be used to verify 
if the system is jammed or not [2o| . 
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FIG. 19: The distribution of forces, P(F/(F)) vM o 

The force distribution for vMD ensemble, 
Pdyn(F/(F)dyn) is shown in Fig. IT8k . From Fig. 
[TSk . we find that the force distribution Pdyn(F / (F) d yn ) 
as a function of different volume fraction 4> no longer 
collapse. At <j> close to <f> c , the average system force F for 
each packing changes dramatically. While at (j> is much 
above <j> c , the fluctuations of the average system force 
F decrease, then the force distribution Pdyn{F / (F) dyn) 
changes continuously. 

We can also calculate the force distribution 



FIG. 20: Microcanonical calculations. The entropy surface 
S(\n((j} — cj> c ),]np). The color bar indicates the value of the 
entropy. The superimposed blue O is (p(<W)dyn from MD 
calculations as in Fig. 1111 The olive arrow line indicates the 
maximization direction of the entropy (— sin 6, cos 9). Fol- 
lowing this direction, the entropy is maximum at the point 
(ln((</>)d yn — 4>c), ln(p)dyn), corroborating the maximum en- 
tropy principle. 



VIII. THERMODYNAMIC ANALYSIS OF THE 
JAMMING TRANSITION 

So far we have considered how the angoricity deter- 
mines the pressure fluctuations in a jammed packing at 
a fixed <p. The role of the compactivity in the jamming 
transition can be analyzed in terms of the entropy which 
is easily calculated in the microcanonical ensemble from 
the density of states. Figure [2D] shows the entropy of the 
system as a function of (p, 0) in phase space: 

S = ln(n(p, <!>)). (33) 

Here O is the number of states which is the unnormalized 
version of g(T,<j)). It is important to note that Fig. [2UJ 
shows the non-equilibrium entropy, in the Edwards sense. 
At the Edwards equilibrium, the entropy is maximum 
respect to changes in <f) and T. We will now see how 
the jammed system verifies the principle of maximum 
entropy. 
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We analyze the entropy surface S(ln(0 — C ), lnp) plot- 
ted versus (ln(0 — c ),lnp) in Fig. [501 When we plot 
superimposed the MD-obtained curve (j?(</>))dyn we see 
that the MD values pass along the maximum of the 
entropy surface constrained by the coupling between 
p and </>, Eq. ((U (such a curve is superimposed to 
the entropy surface in Fig. BO)) . Due to the coupling 
through the contact force law, the maximization of en- 
tropy is not on p or <j) alone but on a combination of 
both. The entropy S reaches a maximum at the point 
5(ln((0)(j yn — 0c), l n (p)d yn ) when we move along the di- 
rection perpendicular to the jamming curve (p(</>))d yn 
(see the maximization direction in Fig. [20]) . This is a 
direct verification of the second-law of thermodynamics: 
the dynamical measures maximize the entropy of the sys- 
tem. 

We can use this result to obtain a relation between an- 
goricity and compactivity and show how a new "jamming 
temperature" Tj and the corresponding jamming "heat" 
capacity Cj can describe the jamming transition. 

From the power-law relation p = T/V oc (0 — 4> c ) a , we 
have: 



Inp = lnpo + a1n((p - (j> c ), 



(34) 



where po is the constant depending on the system. 

Figure [20] indicates that the jammed system always 
remain at the positions of maximal entropy, 



6S = 0, 



(35) 



in the direction (— sin 9, cos 9), perpendicular to the jam- 
ming power-law curve and the slope 



tan 9 — a. 



(36) 



In order to further analyze this result, we plot the en- 
tropy distribution along the direction (— sin (9, cos #) in 
Fig. [2U We see that the entropy of the corresponding 
jammed states remains at the peak of the distributions 
along (— sin 9, cos 9). This is clear when we plot the value 
of (p, <fi) from MD simulations in the plot of S in Fig. |2"T1 
blue dot. Except for volumes very close to jamming, the 
MD coincides with the maximum of S when taken along 
(— sin 9, cos 9). We notice that the maximization is quite 
accurate for large volume fractions. For close to jam- 
ming deviations are seen. We cannot rule out that these 
deviations are finite size effects. The deviations for small 
</> (Fig. [2Tj) remains to be studied. They could be due 
to finite size effects or due to the fact that the value of 
(f> c is different for the MD results and the microcanoni- 
cal ensemble S due to the small size of the system. In 
general, this plot verifies the maximum entropy princi- 
ple in this particular direction. An analogous plot where 
the entropy is shown as a function of but along the 
horizontal direction (or along the vertical direction, F) 
shows that the MD entropy is not maximal along these 
two directions. 

Thus, the maximization of entropy is not on T or V 
alone, but on a combination of both. This means that the 



entropy S , (ln((0)d yn — 4>c), ht(p)dyn) is maximum along 
the direction of (— sin 0, cos 0) and the slope for the en- 
tropy of the jamming power-law curve along this direc- 
tion (— sin^jCos^) is (see Fig. |2"2"]) . that is, 



OS 



dhx((j> - 4> c ) 



sm( 



dS 
d lnp 



cosf 



(37) 




FIG. 21: The non-equilibrium entropy S^lnp, ln(0— <j> c )) along 
the direction (— sin 0, cos 0) for different jamming ensemble 
points. The blue O represents the entropy of the jammed 
system obtained from MD. We see that closely follows the 
maximum of S for all the volume fractions except very close 
to the jamming point where the blue point does not coincide 
with the maximum of S. It remains to be studied if this 
deviation is a finite size effect, or it could be due to a different 
value of (j> c between simulations and microcanonical ensemble. 
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FIG. 22: The representation of the maximization analysis 
SS — along the direction (— sin#,cos#) for one point in 
the jamming power-law curve. Here ci = T and C2 = ((/> — 

0c)(A%/<A 

By the definition of angoricity A = dT/dS and com- 
pactivity X = dV/dS, we have: 



dS 



dS 



P- 



d lnp dp 



,95 

or 



£ = ci 

A A' 



(38) 
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NV g 1 _ c 2 
6 2 X X' 



(39) 



where <f> = NV g /V, c r = T and c 2 = (</» - c )(A%/c?!> 2 ). 
By Eq. J35J and Eq. d35J), we can simplify Eq. (|37|): 



(40) 



The relation between X and A can be obtained then 
(Fig. 



X = -a—A = -a- 
ci 



From Eq. (|4"Tj) we obtain that: X oc 
and near <b r : 



X 



-A. 



(41) 



(42) 



We notice that the compactivity is negative near the 
jamming transition. A negative temperature is a general 
property of systems with bounded energy like spins [H[ : 
the system attains the larger volume (or energy in spins) 
at <fi c when X — > 0~ and not X — > +oo [The bounds 
<j) c < (p < 1 imply that the jamming point at X — > 0~ 
is "hotter" than X — > +oo. At the same time A — » + 
since the pressure vanishes] . 

We conclude that, A and X alone cannot play the 
role of temperature, but a combination of both deter- 
mined by entropy maximization satisfying the coupling 
between stress and strain. Instead, there is an actual 
"jamming temperature" Tj that determines the direction 
(— sin 9, cos 9) in the log — log plot of Fig. [20] along the 
jamming equation of state (see Fig. I22j) . By maximizing 
the entropy along this direction we obtain the "jamming 
temperature" Xj as a function of A and X: 



1 c x 
—— == — — sin 
Tj A 



That is: 



% cos 9 = cos 8(a% - %). (43) 
X A X 



Asin6» 

Cl 



X cos sin ( 



C 2 



a A 



r 

s-y-a 



(44) 



Thus, the temperature vanishes at the jamming transi- 
tion. 

Furthermore, the "jamming energy" corresponding 
to the "jamming temperature" Tj in Eq. (1431) . has the 



relation as below 
d£j = TjdS 
= Tj 



OS 



dln(4> — 4>c 



-dln(</>- 



Tj— dlnp 

o mp 



. X cos 8 . . c 2 . . , , , , , A sin 9 a , , 

( -din - C + T dlnp 

c 2 X c x A 

cos#dln(0 — 4> c ) + sin^dlnp 

(cos 9 + sin 9 tan 8)d ln(</> — <j) c ) 

dln(</> - i^ c ) 



cos 9 



That is, 



and 



dEj = Vo2+TdhM>-0 c ), 



= (Va 2 + l)M0-<Ac). 



(45) 



(46) 



(47) 



By the definition of "heat" capacity, we obtain two 
jamming capacities as the response to changes in A and 
X: 



C T = dT/dA ~ {<t> - 0c)- 1 ~ A- 2 / 5 , 



(48) 



C v = dV/dX ~ (0 - ^e)" 1 ~ lAr 1 ^. 
The jamming capacity Cj can be obtained as: 

r = T — = T -^^l+T — dln ^-^ c 
J J 9Tj J 9 \np dTj 3 8 ln(</> - <j> c ) dTj 

(49) 

Finally, with Eq. ([3"T ]) - (j3"9"|) . the capacity Cj can be cal- 
culated: 

n m /Cl c 2 51np 1 + a 2 c t ding 

aj = _ ^ } otT = T3 -^~Tm' (50) 

Since Tj ~ (4> — </> c ) and p ~ (0 — C ) 15 , we obtain 

Cj^C^-^c)- 1 . (51) 

From Eq. I|48p. the jamming capacities diverge at the 
jamming transition as A — > + and X — > Q~. However, 
this result does not imply that the transition is critical 
since from fluctuation theory of pressure and volume [58[ 
we obtain: 



((Ar) 2 ) = A 2 Cr ~ A i e , 
({AV) 2 )=X 2 C V ~ IAI 1 ' 5 . 



(52) 



Thus, the pressure and volume fluctuations near the jam- 
ming transition do not diverge, but instead vanish when 
A — > + and X — > 0~. From a thermodynamical point 
of view, the transition is not of second order due to the 
lack of critical fluctuations. As a consequence, no diverg- 
ing static correlation length from a correlation function 
can be found at the jamming point. However, other cor- 
relation lengths of dynamic origin may still exist in the 
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response of the jammed system to perturbations, such as 
those imposed by a shear strain or in vibrating modes 
@, Such a dynamic correlation length would not 

appear in a purely thermodynamic static treatment as 
developed here. We note that static anisotropic packings 
can be treated in the present formalism by allowing the 
inverse angoricity to be tensorial (37j . 

The intensive jamming temperature Eq. (144)) gives use 
to a jamming effective energy Ej as the extensive vari- 
able satisfying Xj = dEj/dS and a full jamming capac- 
ity Cj ~ (4> — </> c ) _1 , which also diverges at jamming. 
However, the fluctuations of Ej defined as ((Ai?j) 2 ) = 
TjCj ~ Tj has the same behavior as the fluctuations of 
volume and pressure, vanishing at the jamming transi- 
tion T 3 -> 0+ [A -> 0+ in Eq. (|44)) ]. 



IX. COMPARISON WITH O'HERN ET AL. 




10" 3 



2000 4000 



FIG. 23: Sampling probability of each microstate ft identified 
by its rank k fro low to high. Results are for a system of 30 
particles at <j> = 0.61 and a narrow set of pressures around 0. 



The results so far show a general agreement between 
MD and the ensemble average. These include the maxi- 
mum entropy principle and ergodicity. We now turn to a 
comparison with similar simulations done by O'Hern et 
al. (25l 26J. These studies perform an exhaustive search 
of all configurations in the PEL of frictionless particles 
similarly as in the present paper. However, they find that 
the microstates are not equiprobable, i.e., microstates 
with the same pressure and volume fraction (pressure is 
fixed at zero since only hard sphere states are of inter- 
est) do not have the same probability when sampled by 
a given algorithm. Furthermore, experimental studies of 
equilibration between two systems [28j . suggests that a 
hidden variable is necessary to describe the microstates, 
further supporting the results of [25[ . The applicability of 
the microcanonical ensemble is based on the fact that the 
microstates are defined by (r, 0). Thus, the fact that the 
states are not equiprobable implies that there must be 
an extra variable needed to describe their probabilities. 
Therefore, ergodicity and the maximum entropy princi- 
ple, which are downstream from equiprobability, are not 
supposed to hold, in disagreement with the results shown 
in the present paper. 

To investigate this situation, we repeat the same cal- 
culations as in 25] with our algorithms. We first rule 
out subtleties related to algorithmic dependent results 
in sampling the space of configurations. We use our 30 
particles system and use <f> = 0.61 very close to jamming 
and r = to look for the hard sphere packings. We 
search for the jammed configurations as above. We re- 
call that the sampling of the space of configurations is 
not complete due to the relatively large system size but 
represent a good sampling as discussed above. Ref. [25[ 
uses a different system of 14 particles in 2d for which 
248,900 configurations are found exhaustively sampling 
the phase space (which is estimated to have ~ 371,500 
states). These simulations correspond to a system with 
periodic boundary conditions for which a larger space is 
expected than the close boundary-system of Section III Fl 



However, these differences do not affect the conclusions 
below. 

We start by measuring fk which is the probability to 
find a given microstate k as defined by (2f| : each packing 
can be obtained many times during a search and there- 
fore fk measures the probability for which each packing 
occurs. The main result of [25| is that fk differs by many 
orders of magnitude for states with fixed (T, <f)). Indeed, 
even configurations which are visually very similar can 
be 10 6 more frequent, see Fig. 1 of 25]. 

Figure |2"31 shows fk sorted as a function of k, the rank, 
as in [25]. This plot reproduces the results of [25[ in 
our system. For a fixed pressure and volume there are 
many states with a large difference in their probability. 
The least probable states are 10 -3 less probable than the 
most probable state showing a breakdown of equiproba- 
bility. The question is how to interpret the results of er- 
godicity in the light of the failure of equiprobability and 
whether there is a need for an extra variable to describe 
the microstates. 

We first mention the issue of the small system size. It 
is quite possible that the low probability states will com- 
pletely disappear in the thermodynamic limit and the 
ones remaining are the most probable ones with equal 
probability. Indeed, the flat average assumption is only 
valid in the thermodynamic limit and simply says that 
even if there exists less probable states (10 -3 less proba- 
ble) then they will be irrelevant in the ensemble average, 
thus only the most probable and flat states are impor- 
tant. 

We have done simulations with N =14 particles and 
found that the least probable states are 10 -5 less proba- 
ble than the most probable states. Comparing with the 
factor 10 -3 for N = 30, may indicate that the system 
size may take care of the non-equiprobability problem. 
However, calculations for larger system to fully test this 
assertion are out of the range of current and near future 
computational power. 

Second, we notice that the coordination number is also 
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(b) 

FIG. 24: (a) Sampling probability of each microstate fk as 
a function of the coordination number Zk of each microstate. 
(b) Plot of ln(X flxZfc fk/fJP ax ) versus Z k showing an expo- 
nential decay consistent with the density of states proposed 
in [2<|. 

important to define the jammed states. Figure [M] plots 
the same states as Fig. [23] but as a function of Z^, the 
coordination number of microstate k. The most probable 
states satisfy: 

h{Z k )~e-* z «. (53) 

Furthermore, if we sum up all the states for a given Zk 
and plot log(X]fi x z fc f k ) vs %k we obtain Eq. ([53")) as seen 
in Fig. [24b . This result does not mean that Zk is the 
hidden variable but rather Eq. (|53[) provides the density 
of states proposed in [29| in the thermodynamics calcu- 
lation of the random close packing of spheres. Indeed, 
we have predicted that the density of states g(z) = h z z , 
with h z playing the role of a Planck constant defining the 
minimum size in the volume landscape. According to Eq. 
([53"[). this prediction is satisfied in average with h z = e~ 8 
which is a small number as expected. 

This result indicates that some variability in the prob- 
abilities of the microstates is expected from the fluctua- 
tions in the coordination number of each microstate. In 
Appendix lAl we elaborate an extension of the framework 
of 1 29] to incorporate fluctuations in Z that are neglected 
in [29(. The purpose is to test whether the RCP and jam- 



ming transition are affected by these fluctuations. We 
find that the results are consistent with those found in 
29]. 

We notice that for a fix Zk there are still many 
marginal states with very small probabilities as seen in 
Fig. I2"4"k . If these states do not completely disappear 
in the thermodynamic limit, then they need to be ex- 
plained. We end this discussion by providing a possible 
explanation for the existence of these states. 

The numerical breakdown of equiprobability might be 
related to the fact that the found packings are not indis- 
tinguishable. Indeed, we ignore the rotation and trans- 
lation symmetries of the packing in order to make the 
numerical search possible. However, for the Edwards 
flat hypothesis, these packings should be assumed differ- 
ent. Once we breakdown the rotational symmetry, there 
would be many similar packings. The high degeneracy of 
the high symmetric packings may be responsible for the 
uneven distribution, which would be, in this case, simply 
artificial. 

For instance, consider two packings with 4 particles: 
(a) a square packing with each particle on the corner 
and (b) a triangle with each particle in each corner plus 
one in the center. 

For both packings there are 4! =24 different per- 
mutations, which should be considered as 24 different 
packings, in principle. However, since we can rotate the 
square packing by 90 degree and obtain the same one, 
there are only 24/4 = 6 distinguishable packings. Simi- 
larly, for the triangle, there are 24/3 = 8 distinguishable 
packings. The probability between (a) and (b) is uneven 
(6:8) if we assume that each distinguishable packing is 
equal-probable. Therefore, different symmetries of the 
packings may contribute to the unequal probabilities that 
we measure in the algorithms. 

Therefore, if the Edwards assumption is correct, fk 
should be proportional to Sk, where Sk is the order of 
the symmetry group (point group) of the packing k, since 
there are Sk degenerations (same packing if particles are 
identical). This conjecture needs extra evaluation of the 
symmetry of each packing. For instance, the translation 
invariance is important, and for cubic periodic boundary, 
it is also important to include the symmetry of cubic 
point group C 3h . 

We do not investigate this co njecture but ra ther pro- 
vide the codes and packings in http://jamlab.org to do 
that. Since the 3d case is complicated, one might try the 
2d system first to easily visualize different packings. A 
simple question is: given two packings with different fre- 
quencies, how do they look like (25[? Would be the high 
symmetric one visited more, or inversely? 



X. CONCLUSION 

We have demonstrated that the concept of " thermal- 
ization " at a compactivity and angoricity in jammed sys- 
tems is reasonable by the direct test of ergodicity. The 
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numerical results indicate that the full canonical ensem- 
ble of pressure and volume describes the observables near 
the jamming transition quite well. From a static thermo- 
dynamic viewpoint, the jamming phase transition does 
not present critical fluctuations characteristic of second- 
order transitions since the fluctuations of several observ- 
ables vanish approaching jamming. The lack of critical 
fluctuations is respect to the angoricity and compactivity 
in the jammed phase <fi — > which does not preclude 
the existence of critical fluctuations when accounting for 
the full range of fluctuations in the liquid to jammed 
transition below (j) c . Thus, a critical diverging length 
scale might still appear as (j> — > d >~ J60l| . which has been 
recently observed by experiment [6l|. 

In conclusion, our results suggest an ensemble treat- 
ment of the jamming transition. One possible analytical 
route to use this formalism would be to incorporate the 
coupling between volume and coordination number at the 
particle level found in [2Sj, |62[ together with similar de- 
pendence for the stress to solve the partition function. 
This treatment would allow analytical solutions for the 
observables with the goal of characterizing the scaling 
laws near the jamming transition. 
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Appendix A: Microstates and Fluctuations in 
coordination number 



Here, we develop a Z-ensemble for hard spheres in 
the limit of zero angoricity. In the main test we found 
that fluctuations in Z may account for certain variabil- 
ity in the probability of microstates. Here we investigate 
whether this variability affect the existence of RCP and 
the jamming point. We develop a partition function in 
Edwards ensemble to study the dependence of RCP on 
this type of fluctuations. 

The partition function is 



Nz 



JJ e -(*/**+/s«/*i)efe. j 



(Al) 



where z m m = Z and z max = 6, f3 = 1/X, and k = 2y3. 
We follow the notation and concepts from [2!| |62T - [64| . 
We define x = {J2i z i)/N, thus: 



where 
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P(x)dx, 
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1 (A3) 

where z* — 1/8 according to Fig. I2~ib . We consider the 
inverse Fourier transform of P x (f): 
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Since 



Then: 
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where K n (a) is the modified Bcssel function of the second 
kind. By taking the coupling constant 
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a = 2B 1 / 2 , 
z = B 1/2 z*x. 
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Thus. 
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Now, we expand 



ln(I + x) = x — -\-x 2 + ]-x 3 
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is just a Gaussian distribution with the mean 
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Thus, by using the saddle point approximation, we ob- 
tain the free energy density /: 
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We also obtain the energy density, or volume density 
in the context of Edwards: 
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k rfB " 25 K {2B 1 / 2 ) 
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wllCre Z m ax — Zmaxfz •> Z m in = Zmin/z &nd 



because 



L(B) = B K^B^-y 
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Zmax) 
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Thus. 



(3w = 



L(B) 



+ 



1 (L(B){L{B) 

2 V B + L(B) - L(B) 2 

1 / L(B)(L(B) - Z min ) 

2 V B + L(B) - L(B) 2 



{B - Z max L(B))(L(B) 

B + L(B) — L{B) 2 

(B - Z mm L(B))(L(B) — Z min ) 
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Q(L(B) - Z max ) 
Q(Z min - L(B)), 



(A22) 



and the entropy density: 

s = P(w-f). (A23) 

There are two phase transitions at L(B) = Z m i n and 
L{B) = Z max . For the jammed phase Z m i n < L(B) < 
Zmax, we have f3w = L(B) — 1/2. If z* is a small value, 
z* = 1/8 from Fig. [Ml then B is relatively large. Thus, 
L{B) w B 1 ' 2 and w rnax w L(B)//3 = k/z*B- 1 / 2 = 
n/(z*Z min ) = n/z min . Similarly, w min w n/z rnax , which 
is consistent with the boundaries of the phase diagram 
obtained in [2!|. Furthermore, / sa 2B 1 / 2 and s w 



s - B 1/2 , where s = Z mM . Or, s = (z max - n/w)/z*. 
Thus, we have verified that the inclusion of fluctuations 
in the coordination number does not change the shape of 
the jamming phase diagram obtained in [29, 63]. These 
fluctuations may affect the probability of the microstates 
according to the density of states proposed in [2^]. A 
further application of this generalized Z-ensemble is de- 
veloped in 65] to calculate the probability of coordina- 
tion numbers in packings, with good agreement with the 
numerical results for different packings in the phase dia- 
gram. 



